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In this paper, we present a theoretical study of a Bose-Einstein condensate of interacting 
bosons in a quartic trap in one, two, and three dimensions. Using Thomas-Fermi approx- 
imation, suitably complemented by numerical solutions of the Gross-Pitaevskii equation, 
we study the ground sate condensate density profiles, the chemical potential, the effects of 
cross-terms in the quartic potential, temporal evolution of various energy components of 
the condensate, and width oscillations of the condensate. Results obtained arc compared 
with corresponding results for a bosc condensate in a harmonic confinement. 
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1. Introduction 

Since the first experimental observation of Bose-Einstein (BE) condensatiorfll^IS] 
in bose atom clouds in a harmonic trap, both experimental and theoretical studies 
of this phenomenon and variou s p roperties of the condensate has grown up to a 
fast expanding field of research SB 6 7 8 9 10l_ \Yhile studies of Bose condensation 
in a harmonic potential is of great interest, investigations of the BE condensation 
of lattice-bosons and free-bosons in anharmonic potentials, in p articular in a q uar- 
tic potential, have received considerable attenti on in 

recent timeAllMMiel 

see 

also the note in Ref. . All the previous workJUIEIllIIl] ^ except twJlSIlS]^ deals 
with lattice-bosons in anharmonic confinements. The previous wor ks on f ree-bosons 
in power law traps and quartic traps in three and lower dimension j-^^ l -^^ -^ dealt with 
the calculations of Bose condensation temperature for non-interacting bosons in 
those environments. In this paper, we consider interacting bosons in a quartic po- 
tential. The study of interacting bosons in a quartic potential may have ir npo rtant 
implications for fields ranging from chaotic dynamics to high energy physical^. Fur- 
ther, it would be desirable to have minimum inhomogeneity in the density profiles 
for clear observation of several phases and tr ansitions possible for interacting bo- 
son, fermion, and boson-fermion mo 

deiPini. 

While interactions in a many-boson 
system are expected to broaden the real space distribution of bosons in the confin- 
ing region, it is not clear how far, from the center, a homogeneous distribution of 
bosons will exist. We will show, among other results, that the interaction leads to 
a more or less homogeneous boson distribution in a substantial region (compared 
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to the non interacting bosons case) in and around the trap center. It would also 
be of contemporary interest to investigate the dynamics of the condensate in the 
specific quartic potential environment and compare it with the usually employed 
quadratic environment. It may be expected that the experiments would first involve 
investigations of interacting bosons in a quartic potential before an additional lat- 
tice potential is applied. Furthermore, experimental studies of interacting bosons 
in the simplest anharmonic environment may be of considerable interest in its own 
right. Within this general setting, it would be useful to make a detailed theoretical 
study of interacting bosons in a quartic potential. It is the purpose of this paper to 
present some results of such an investigation using a combination of analytical and 
numerical methods. 

This paper is organized in the following three sections. In Sec [21 the conden- 
sate properties are analyzed through a stud y of th e Gross-Pitaevskii equatiorP^IIS 
employing the Thomas- Fermi approximatiorPOEl] jj^ [3] present a compara- 
tive study of real space distributions of bosons in condensates of interacting bosons 
in quartic and harmonic confinements, by numerically solving the Gross-Pitaevskii 
equation. The time evolution of some properties (various components of energy and 
width oscillations of the condensate) of the condensate is given in Sec IH and the 
conclusions are presented in Sec [3 



2. BEC in a quartic trap: Thomas-Fermi approximation 

The ground state properties of a Bosc-Einstcin condensate (BEC) is well described 
by a macroscopic wave- functio n ^jfr, t ) , the time evolution of which is governed by 
the Gross-Pitaeveskii equation SSHSl (GPE). Incorporating the trap potential as 
well as a repulsive interaction between the atoms forming the condensate, the GPE 
in three dimensions (3d) is given by 



= -^^VV(r, t) + «(r)^(r, t) + iVf/o| V(r, t)P^(r, t) , (1) 



dt 2m 

where r = {x,y,z)'^, m is the atomic mass, t;(r) = Q!(r.r)^ is the quartic trapping 
potential, N is the total number of atoms forming the condensate, Uo{= Airfi^a/m) 
is the strength of interaction between atoms in the condensed state, and a (> 0) is 
the s-wave scattering length. 

Before applying the Thomas-Fermi approximation, we non-dimensionalize Eq. 
([T]) through a set of linear transformations: i = ujt, f = r/l, ■ip{r) = /■^/^'0(r). After 
dropping the wiggles on the symbols, we obtain 

'^^Ir^ = -^V2V(r, t) + 7 (r.r) V(r, t) + X N \^{r, t)pV>(r, t) , (2) 

where 
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In the absence of the non-linear interaction term, the GPE is reduced to a 
Schrodinger equation (with a quartic potential) whose approxima te ei genvalues can 
be obtained through perturbation theory or WKB approximation Expressing a 
in units of m?LO^/h, we perform dimensional reductions to obtain correspond- 
ing GPEs in lower dimensions. In general, the GPEs for d = 1, 2, and 3 can be 
summarized in a single equation as 

'^^^ - -^V2^(r,t) + (r.r)>(r,i) +ft,|^(r,OlV(r,t) , (4) 

where 

kd = adXdN, Xd = -jjz^- (5) 

Here a2 and ai are scaling factors introduced to write GPE in 2d and Id, respec- 
tively. In order to find a stationary solution of Eq. , we do a separation of variables 
tjj{r,t) = fpir) X exp[—i{^/{huj))t], where /i is the chemical potential. Starting from 
Eq. dU, we obtain 

- ^VVW + ^VW + fedlV-WlVW = -^^{r) . (6) 

In the preceding equation r** = x'^ + + z^. For simplicity we have neglected the 
cross terms from the potential (see, however. Sec. III-C where these terms are shown 
to have only minor effects on the properties we study) . To analytically solve for the 
ground state of the stationar y state GP E, as given in Eq. ([6]), we make use of the 
Thomas-Fermi approximationSEOElEl] -p^j, sufficiently large number of atoms, the 
kinetic energy is very small compared to the potential and the interaction energies. 
So, a good approximation to the ground state is obtained by solving the GPE 
without the kinetic energy term. On using this approximation, the GPE becomes 

r4V'(r) + fcd|^(r)|V(r) = -^V(r), (7) 
The relevant solutions of the preceding cubic equation are 



Since ipii^) vanishes for the Thomas- Fermi radius 7'xf is obtained as 

in which the chemical potential /ixF is be determined from the normalization con- 
dition 

Cd |^(r)|2r^-idr = l, (10) 

Jo 

where Cd = 2, 27r, and An for d = 1, 2, and 3, respectively. Then, the chemical 
potential is found to be 



4:Cd ' 



(11) 
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So, for the 3d case, the chemical potential varies as Nr for the quart ic trap . In 
comparison, we note that it varies as TVt for bosons in a harmonic trajPSEll (j^j. 
which v{r) oc r.r). It may be noted that the N dependence of the chemical potential 
becomes more dominant as the dimensionality decreases. 

The scaling factors {a2 and ai) were introduced in 2d and Id GPE to obtain 
GPEs in 2d and Id. They also ensure that the size of the condensate remains the 
same in all the three dimensions. To find a2 and ai, we c alculate ^tf in 2d and Id 
and equate it to the chemical potential obtained in 3cPS. Using this procedure, we 
obtain 



and 



as = 1 , 



(12) 
(13) 

(14) 



3. Spatial distribution of bosons in the condensate 
3.1. Condensate density profiles in Id 

In this section, we present results on the spatial distribution of bosons in the con- 
densate in Id by numerically solving the G PE. We solve the GPE by using the finite 
difference Crank-Nicholson (FDCN) methocPS125| j^^ ^^lis method for the GPE given 

by 

I — — — = Hip[x,t), (15) 

the solution is advanced in small time steps St as given below 

1 - i§lH 

i^{x,t + 6t) = .l i'{x,t). (16) 

i + ^~2^ 

The ground state solution is found by the imaginary time propagation method^. 
In the imaginary time propagation method, we replace 5t by ~i5t and propagate 
a initial trial wave function using the above mentioned propagation scheme. After 
suitable normalization at each time step the wave-function converges to the ground 
state. We have taken the space step as 5x = 0.1 and the time step as 5t = 0.001. 
For illustrative purpose, we have used a set of parameters for ^^Rb: m — 1.44 x 
10-25 iCg, a = 5.1 X 10^^ to, v {= uj/2tx) = 2AHz. 

In Fig. 1, we show the condensate density profiles for systems with 10 to 10^ 
bosons in a quartic trap. Further, this figure also contains results of the TF approx- 
imation for the same numbers of bosons in a quartic trap. Furthermore, we have 
also shown the results for bosons in a harmonic potential. We find that the TF 
approximation results are in quantitative agreement with the numerical solutions 
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of GPE when the number of bosons is large. On increasing the number of atoms, 
the width of the g round state increases, which is similar to the case of bosons in a 
harmonic 

potentiaPSm^ 

except that the spread of condensate is larger for a har- 
monic potential. We also notice that the ground state density profiles for the quartic 
trap becomes relatively flat with increasing number of bosons in the quartic trap. 
Next, we fix the number of atoms to 10^ in a quartic trap and investigate the boson 
distribution in the ground state for different interaction strengths. The results are 
displayed in the left panel of Fig. 2. We find that the ground state boson distri- 
bution is more spread out with increasing interaction strength. It may be argued 
that the redistribution of bosons from the center to other parts of the condensate 
is due to the peak in the non-interacting distribution in the central region since the 
U = wave-function is peaked at the center and is symmetric around the center. 
But. for the interacting bosons case, the interaction term is to be calculated in a 
self-consistent manner including the effect of the interaction. Hence, it is not obvious 
that interacting bosons distribution will have the profile shown in Fig. 2. It is the 
result of a combined effect of interaction and self-consistent changes in the wave- 
function. Thus, we see that the interacting bosons, compared to the non-interacting 
bosons, have a density profile which is more or less uniform in the central region 
in the case of a quartic confining potential. In Fig. 2, we have also exhibited the 
ground state chemical potential for different values of the interaction parameter fci . 
The chemical potential for a condensate in a quartic trap is larger than that in 
harmonic trap, the difference becoming more significant for higher values of ki. 



3.2. Condensate density profiles in 2d and 3d 

We have solved the GPE in higher dimensions by using the split-step FDCN 
methocP^. In this method, in 3d, the Hamiltonian H is split into different non- 
derivative and derivative parts Hi, H2, i?3, and H4, where 

i/i =i.(r)+fc3|V'(r,t)l', (17) 



1 92 1 52 Q2 

The time evolution is performed in the following steps. Let '0" be the wave-function 
at time i„ . This wave- function is advanced first over a time step 6t at tn to produce 
an intermediate solution from -0" via 

Next the time propagation is performed via the following semi-implicit Crank- 
Nicholson scheme: 

r^'" - TT^^"""'^'' (20) 
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,n+3/4 




,n+2/4 



(21) 



and 



1 



1 + i 




,n+3/4 



(22) 



An analogous set of equations hold in 2d. In the calculations, we have used space 
steps of 0.1 (for Sx and Sy in 2d and Sx, Sy, and Sz in 3d) and the time step 
St — 0.001. To calculate the ground state, again we use imaginary time propagation 
method by replacing St with —iSt. We start with an initial trial wave-function which 
is propagated via the above mentioned scheme. After suitable normalization, at each 
time step, the wave-function converges to the ground state. 

In Fig. 3, the ground state density profiles of the condensate are shown for 
different numbers of atoms ranging from lO'' to 10^. On increasing the number of 
atoms, the width of the ground state increases, which is similar to that in Id. On 
comparing the ground states for harmonic and quartic traps as shown in Fig. 4, we 
found that the ground state density profiles are flat in the central region for the 
quartic trap compared to the harmonic trap. Similar results are obtained for the 3d 
case as well. Clearly, the quartic confinement is more advantageous to the quadratic 
one if one is interested in investigations of the properties of a bose condensate with 
a minimal influence of the confining potential. We may note here that, in an earlier 
wor we have shown that the finite temperature properties of lattice bosons are 
closer to the pure lattice case if one uses quartic rather than a quadratic confining 
potential. The effect of changing interaction strength on the ground state density 
profiles in a quartic trap are displayed in Figs. 5 and 6 for 2d and 3d, respectively. It 
is clear that the bosons gain energy by spreading into larger area in the traps with 
increasing strength of the interaction. In Fig. 7, the ground state chemical potential 
is plotted for different values of the interaction parameter kd for d = 2 and 3. The 
chemical potentials of the condensate in a quartic traps are higher than those in 
harmonic traps, the difference becoming more significant at higher values of fc^. 

3.3. Effects of cross terms 

In the preceding section, we have presented results of our calculations neglecting 
the cross terms in the full quartic potential {V{r) ~ a{r.r) ). In this section, we 
investigate the effects of the previously neglected cross terms. The contribution of 
the cross terms is in cross-coupling potentials across different dimensions. One may 
expect that the inclusion of these cross-terms may lead to complex inhomogencous 
density profiles. But, on performing the calculations for the quartic trap, we find 
that the effects of cross-terms are rather minor. They lead to small increases in 
the chemical potential and the number of bosons in the central region of the trap 
as shown in Fig. 8 for 2d. We do not find any significant inhomogeneities in the 
ground state condensate profiles generated by the cross-terms. We have observed 
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(figures not shown) similar effects in 3d as well. It may be noted here that cross 
terms lead to a 10% increase of the bose condensation temperature (Tc) in 2d, and 
a 20% increase in 3cff^. 



4. Time evolution of the condensate 

In this section, we study the time evolution of the condensate after removal of 
the trapping potential and width oscillations of the condensate. To theoretically 
study the expansion of condensate, we time evolve the GPE after removing the 
potential term, taking the stationary ground state as the initial state function at 
t = 0. To time evolve the GPE we use the FDCN scheme given by Egs. [TMTHl in Id 
and split-step FDCN scheme as given in Eqs. [I7ll22l in higher dimensions. Here the 
propagation is done in real time with time-step 5t = 0.0001 and space steps of 0.1 
(for Sx, Sy, and 6z). The different energy components we calculate are given bjS^^Hl 

^kin(t) = -l f r{r,t)Vmr,t)dr, 



2 „ 

Epot{t) = j r{r,t)vir)i;{r,t)dr, 

Eintit) = \j r{r.t)kd\^{r,t)\^^{r,t)dr. (23) 
And, the the chemical potential is given 

M(i)- / Q|VV'(r,t)P+«(r)|V;(r,t)|2-f-fc,|V;(r,t)|4^dr (24) 

The variation of various energy components and the chemical potential, for 
bosons in quartic and harmonic traps, with time is shown in Fig. 9 for Id and 
2d. Similar results (not shown) are observed in 3d. As shown in these figures, the 
chemical potential and the interaction energy decreases on time evolution, but there 
is an increase in kinetic energy term. Most of the energy before the expansion {t < 0) 
is contained in the mean-field and the potential energy of the bosons. At i = 0, the 
confining potential is removed and the corresponding energy term from then on 
is zero. During the expansion, the chemical potential decreases due to decrease in 
the mean field. Further, the increase in kinetic energy during the expansion is due 
to transfer from the mean-field. We have also shown in Fig. 9, the release energy 
(E'rci = -Ekin + Eint) which is a a time independent quantity. 

The time evolution of the width of the condensate is displayed in Fig. 10 for Id, 
2d, and 3d. The condensate widths along three axes x, y, and z are defined as Ax = 
a/< {x~ < X >)2 >, Ay ^ y/< (y- < y >)2 >, and Az (z- < z >)2 >, 

where brackets <> represent the expectation value at a particular instant during 
the time evolution. We assume that at i = 0, the condensate is in a non-interacting 
ground state. Once the interaction is switched on at t = 0, it produces oscillations 
in the condensate. The Fig. 10 show that the frequencies of oscillations are greater 
in the quartic potential case than in the harmonic casJ^. For isotropic quartic 
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potential, amplitude and frequency of oscillations in x,y, and z directions are the 
same. 

5. Conclusions 

In this paper, we studied the static and dynamic properties of a BE condensate in 
a quartic trap in one, two, and three dimensions. The analysis was done using both 
analytical and complementary numerical techniques. The analytical solutions were 
obtained using the Thomas-Fermi approximation. The solutions obtained were then 
compared with the results obtained by numerically solving the GPE. We compared 
the ground state density profiles of the condensate in a quartic trap for different 
strengths of interaction by changing the scattering length. It was found that the 
Thomas-Fermi approximation results were in remarkably good agreement with the 
results obtained from the GPE for a large number of atoms in the condensate in 
a quartic confining potential. The results were also compared with corresponding 
results for a bose condensate in a harmonic confining potential. We find that the 
condensate density profiles are relatively flat in the central region of the quartic po- 
tential compared to the case of a condensate in a harmonic potential. An implication 
of this result is that for investigating properties of a bose condensate with minimal 
influence of the overall confining potential, the quartic potential is more preferable 
compared to a quadratic potential. We also find that increasing the number of atoms 
in the condensate or increasing the interactions spreads the condensate in the con- 
fining region. Our calculations of the chemical potential in the ground state show 
that it is larger for a condensate in a quartic trap compared to that in a harmonic 
trap. The preceding conclusions were arrived at without including the cross-terms in 
a general quartic potential. On including the cross-terms, we find that their effects 
lead to increases in the chemical potentials and condensate density in the central 
region of the trap. These changes were found to be rather minor. From our results 
for the time evolution of the various energy components of an expanding conden- 
sate (after removing the confining potentials), we find that the chemical potential 
and the potential energy decrease during the expansion while the kinetic energy 
increases. It was also noticed that changes in the case of quartic potential are larger 
than the case of a harmonic potential. The condensate width oscillations produced 
by switching on of the interaction, starting with a condensate of non-interacting 
bosons, were also studied. It was found that the width oscillations of a condensate 
in a quartic trap have higher frequency than that in a harmonic trap. 
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Fig. 1. The density profiles of the condensate ground state in Id for various values of N: 10 (top 
left panel), 10^ (top right panel), 10* (bottom left panel), 10^ (bottom right panel). The results 
are for; the quartic trap (solid lines), TF results for the quartic trap (dotted lines), and for the 
harmonic trap (dash-dot lines). In this and other figures x is in units of I. The value of Ai is 
1.41 X 10"" m2. 
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Fig. 2. Left panel: The density profiles of the condensate ground state of IC* bosons in a Id 
quartic potential for different values of the interaction parameter fci . The curves are for (from top 
to bottom): fci = 0, 22.52, 30.07, 49.34, 80.95. Right Panel: Ground state chemical potential fi 
for various values of interaction ki in Id for bosons in a quartic potential (solid line) and in a 
harmonic potential (dotted line). The /i is in units of hui and x is in units of I. 




Fig. 3. The density profiles of the condensate ground state in 2d quartic trap for: A'^ = 10^ (top 
left), N = 10" (top right), N = 10^ (bottom left), and TV = lO" (bottom right). The value of A2 
is 6.408 X 10^* m. The x and y arc in units of I. 
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Fig. 4. The density profiles of the condensate ground state in 2d for N = 10* for the quartic 
potential (top) and for the harmonic potential (bottom). The value of A2 is 6.408 X 10^* m. The 
X and y arc in units of I. 




Fig. 5. The density profiles of the condensate ground state in 2d for the quartic potential for 
different values of the interaction parameter k2 ■ The figures are for N = 10^ and for: k2 = 232.25 
(top left), 128.21 (top right), and 70.78 (bottom left), 50.00 (bottom right). The x and y arc in 
units of I. 
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Fig. 6. The density profiles of the condensate ground state in 3d for a quartic potential for 
different values of the interaction parameter fcs. These figures are for A'^ = 10* and for: k-j = 
581.79 (top), 145.44 (bottom). The values of {xo,yo) are: (0.1,0.1) (top), (0.1,1.5) (middle), and 
(0.1,2.0) (bottom). The x, y, and z are in units of I. 




Fig. 7. A comparison of chemical potentials with interaction strengths for bosons in 2d (left 
panel) and in 3d (right panel) quartic (solid lines) and harmonic potentials (dotted lines). The fj, 
is in units of fiu!. 
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Fig. 8. Left panel: The variation of ground state chemical potential with interaction strength 
A;2, in a 2d quartic potential, including cross-terms (dotted line) and without cross-terms (solid 
line). Right panel: The density profiles of the condensate ground state in a 2d quartic potential for 
quartic potential with cross terms (solid line) and without cross terms (dotted line). The results 
shown are for A*^ = 10* and A2 = 6.408 X lO^^m. The ^ is in units of hui and the x and y are in 
units of I. 




Time Time 



Fig. 9. The evolution of chemical potential (/i) and energy components with time in Id (left 
panels) and 2d (right panels). The chemical potential (solid lines), the release energy (double-dot 
lines), the interaction energy (dash-dot lines), the kinetic energy (dotted lines), the potential energy 
(dashed lines). The top panels are for quartic potential and the bottom panels are for a harmonic 
potential. The results shown are for TV = 10*, Ai = 1.41 X lO'^^m^, and A2 = 6.408 X 10~^m. 
All the energies and fi arc in units of huj and the time is in units of w. 
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Fig. 10. Condensate widths (Ax, Ay, Az) as a function of time in: Id (top left), 2d (top right), 
and 3d (bottom). The values of fci, k2, and fca are 2.0. The solid lines are for the quartic potential 
and the dotted lines for harmonic potential. The time is in units of w. 



